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Summary 

We  report  here  on  the  work  performed  during  the  first  year  (october  2009  - 
October  2010)  of  the  contract  FA  8655-10-C-4002  on  Multiscale  problems  in 
materials  science:  a  mathematical  approach  to  the  role  of  uncertainty. 

The  bottom  line  of  our  work  is  to  develop  affordable  numerical  methods  in 
the  context  of  stochastic  homogenization.  Many  partial  differential  equations 
of  materials  science  indeed  involve  highly  oscillatory  coefficients  and  small 
length-scales.  Homogenization  theory  is  concerned  with  the  derivation  of  av¬ 
eraged  equations  from  the  original  oscillatory  equations,  and  their  treatment 
by  adequate  numerical  approaches.  Stationary  ergodic  random  problems 
(and  the  associated  stochastic  homogenization  theory)  are  one  instance  for 
modelling  uncertainty  in  continuous  media.  The  theoretical  aspects  of  these 
problems  are  now  well-understood,  at  least  for  a  large  variety  of  situations. 
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On  the  other  hand,  the  numerical  aspects  have  received  less  attention  from 
the  mathematics  community.  Standard  methods  available  in  the  literature 
often  lead  to  very,  and  sometimes  prohibitively,  costly  computations. 

In  this  report,  we  first  review  an  approach  popular  in  particular  in  the 
computational  mechanics  community,  which  is  to  try  and  obtain  bounds  on 
the  homogenized  matrix,  rather  than  computing  it.  Only  computations  of 
moderate  difficulty  are  then  required.  However,  we  will  show  that,  not  un¬ 
expectedly,  this  method  has  strong  limitations. 

We  will  next  introduce  a  class  of  materials  of  significant  practical  rele¬ 
vance,  that  of  random  materials  where  the  amount  of  randomness  is  small. 
They  can  be  considered  as  stochastic  perturbations  of  deterministic  materials , 
in  a  sense  made  precise  below.  We  will  adapt  to  such  a  case  the  well-known 
Multiscale  Finite  Element  Method  (MsFEM),  and  design  a  method  which  is 
much  more  affordable  than,  and  as  accurate  as,  the  original  method. 

The  works  described  below  have  been  performed  by  Claude  Le  Bris  (PI), 
Frederic  Legoll  (Co-PI),  and  Florian  Thomines  (first  year  Ph.D.  student). 

1  Introduction 

Many  partial  differential  equations  of  materials  science  involve  highly  oscilla¬ 
tory  coefficients  and  small  length-scales.  Homogenization  theory  is  concerned 
with  the  derivation  of  averaged  equations  from  the  original  oscillatory  equa¬ 
tions,  and  their  treatment  by  adequate  numerical  approaches.  Stationary 
ergodic  random  problems  are  one  of  the  most  famous  instances  of  mathe¬ 
matical  uncertainty  of  continuous  media.  However,  the  elaborate  tools  and 
techniques  of  (i)  mathematical  probability,  stochastic  analysis,  and  (ii)  nu¬ 
merical  analysis  and  large-scale  computing  have  not  yet  permitted  practical 
computations.  These  are  most  often  accomplished  otherwise  by  the  engineer¬ 
ing  community,  using  more  traditional  approaches.  Despite  definite  achieve¬ 
ments  by  leading  experts,  numerical  analysis  of  stochastic,  and  more  gener¬ 
ally  speaking  non  periodic,  homogenization  problems  remains  in  its  infancy. 

The  purpose  of  this  report  is  to  present  the  recent  progress  we  have  made 
during  last  year  on  this  topic,  with  the  aim  to  make  numerical  random  ho¬ 
mogenization  more  practical.  Because  we  cannot  embrace  all  difficulties  at 
once,  the  case  under  consideration  here  is  a  simple,  linear,  scalar  second  or¬ 
der  elliptic  partial  differential  equation  in  divergence  form,  for  which  a  sound 
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theoretical  groundwork  exists.  We  focus  here  on  different  practical  compu¬ 
tational  approaches. 

This  report  begins,  in  Section  2,  with  a  brief  introduction  to  stochastic 
homogenization.  There  is  of  course  no  novelty  in  such  an  introduction,  the 
only  purpose  of  which  is  the  consistency  of  this  report  and  the  convenience  of 
the  reader  not  familiar  with  the  theory.  We  will  recall  there  why  stochastic 
homogenization  often  leads  to  extremely  expensive  computations. 

In  Section  3,  we  describe  a  classical  approach  from  the  applied  commu¬ 
nities,  which  is  to  try  and  obtain  bounds  on  the  homogenized  matrix,  rather 
than  computing  it.  The  computational  gain  is  evident.  We  will  report  on 
some  numerical  experiments.  Such  experiments  are  likely  to  not  be  new.  But 
they  at  least  show,  quantitatively  and  qualitatively,  the  strong  limitations  of 
such  an  approach. 

As  pointed  out  above,  random  homogenization  for  general  stochastic  ma¬ 
terials  is  very  costly.  Yet,  it  turns  out  that  it  is  possible  to  identify  classes  of 
materials  of  significant  practical  relevance,  where  stochastic  homogenization 
theory  and  practice  can  be  reduced  to  more  affordable,  less  computationally 
demanding  problems.  These  materials  are  neither  periodic  (because  such  an 
oversimplifying  assumption  is  rarely  met  in  practice),  nor  fully  stochastic. 
They  can  be  considered  as  an  intermediate  case,  that  of  stochastic  perturba¬ 
tions  of  deterministic  (possibly  periodic)  materials.  The  case  when  the  tensor 
describing  the  properties  of  the  material  is  the  sum  of  a  periodic  term  and 
a  small  random  term  is  an  instance  of  such  an  approach.  In  Section  4,  we 
show  that  we  can  adapt  to  that  particular  setting  the  well-known  Multiscale 
Finite  Element  Method  (MsFEM),  which  is  designed  to  directly  address  the 
highly  oscillating  elliptic  problem,  rather  than  studying  the  limit  problem 
when  the  typical  small  lengthscale  goes  to  0.  This  method  has  been  initially 
proposed  for  deterministic  problems  [24,  21,  22,  14],  and  has  been  recently 
adapted  to  the  stochastic  setting  [18].  It  then  leads  to  extremely  intensive 
computations.  We  show  in  the  sequel  that,  if  the  problem  is  only  weakly 
stochastic,  then  it  is  possible  to  design  a  method  as  accurate  as  the  original 
MsFEM,  with  a  much  smaller  computational  cost.  As  we  explain  below,  this 
method  is  accurate  provided  the  stochastic  perturbation  is  indeed  small. 

We  collect  in  Section  5  some  conclusions  about  the  work  performed  so 
far,  and  future  directions  for  the  next  two  years  of  contract. 
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2  Basics  of  stochastic  homogenization 

[Detailed,  presentation  can  be  read  in  [1].] 

Stochastic  homogenization  is  best  understood  in  the  light  of  the  easiest 
context  of  homogenization:  periodic  homogenization.  This  is  the  reason  why 
we  begin  with  Section  2.1  laying  some  groundwork  in  the  periodic  context, 
before  turning  to  stochastic  homogenization  per  se  in  Section  2.2. 

We  refer  to,  e.g .,  the  monographs  [15,  19,  25]  for  more  details  on  homoge¬ 
nization  theory,  and  to  the  review  article  [1]  that  we  wrote,  addressing  some 
computational  challenges  in  numerical  stochastic  homogenization.  A  super 
elementary  introduction  is  contained  in  [13]. 

In  this  section,  we  present  classical  results  of  the  literature.  The  reader 
familiar  with  stochastic  homogenization  can  proceed  directly  to  our  contri¬ 
butions,  detailed  in  Sections  3  and  4. 


2.1  Periodic  homogenization 

For  consistency,  we  recall  here  some  basic  ingredients  of  elliptic  homogeniza¬ 
tion  theory  in  the  periodic  setting.  We  consider,  in  a  regular  bounded  domain 
V  in  Rd,  the  problem 


-div 


±per 


X 


Vue 


=  /  in  V, 


ue  =  0  on  dD, 


(1) 


where  the  matrix  Aper  is  symmetric  definite  positive  and  Zd-periodic.  We 
manipulate  for  simplicity  symmetric  matrices,  but  the  discussion  carries  over 
to  non  symmetric  matrices  up  to  slight  modifications. 

The  microscopic  problem  associated  to  (1),  called  the  corrector  problem 
in  the  terminology  of  homogenization  theory,  reads,  for  p  fixed  in  W1, 

|  -div  (Aper(y)  (p  +  Vwp))  =  0  in  Md, 

1  wp  is  Zd-periodic. 

It  has  a  unique  solution  up  to  the  addition  of  a  constant.  Then,  the  homog¬ 
enized  coefficients  read 


+  Vwei(y))T  Aper(y)ejdy, 


(3) 
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where  Q  is  the  unit  cube,  and  where  wei  denotes  the  solution  to  (2)  for  p  =  et, 
with  et  the  canonical  vectors  of  Wl.  The  main  result  of  periodic  homogeniza¬ 
tion  theory  is  that,  as  e  goes  to  zero,  the  solution  ue  to  (1)  converges  to  u* 
solution  to 


— div  [A*S7u*]  =  f  in  V, 
u*  =  0  on  dT>. 


(4) 


Several  other  convergences  on  various  products  involving  Aper  y— J  and  ue 
also  hold.  All  this  is  well  documented. 


The  practical  interest  of  the  approach  is  evident.  No  small  scale  e  is 
present  in  the  homogenized  problem  (4).  At  the  price  of  only  computing  d 
periodic  problems  (2)  (as  many  problems  as  dimensions  in  the  ambient  space) 
the  solution  to  problem  (1)  can  be  efficiently  approached  for  z  small.  A 
direct  attack  of  problem  (1)  would  require  taking  a  meshsize  smaller  than  e. 
The  difficulty  has  been  circumvented.  Of  course,  many  improvements  and 
alternatives  exist  in  the  literature. 


2.2  Stochastic  homogenization 

The  mathematical  setting  of  stochastic  homogenization  is  more  involved  than 
that  of  the  periodic  case. 

We  put  ourselves  in  the  usual  probability  theoretic  setting  for  stationary 
ergodic  homogenization ,  with  the  exception  that  our  notion  of  stationarity 
is  discrete.  It  intuitively  means  the  following.  Pick  two  points  x  and  y  ^ 
x  at  the  microscale  in  the  material  and  assume  y  =  x  +  k  with  k  £  Ul. 
The  particular  local  environment  seen  from  x  (that  is,  the  microstructure 
present  at  x)  is  generically  different  from  what  is  seen  from  y  (that  is,  the 
microstructure  present  at  y).  However,  the  average  local  environment  in  x 
is  assumed  to  be  identical  to  that  in  y  (considering  the  various  realizations 
of  the  random  material).  In  mathematical  terms,  the  law  of  microstructures 
is  the  same.  This  is  stationarity.  On  the  other  hand,  ergodicity  means  that 
considering  all  the  points  in  the  material  amounts  to  fixing  a  point  x  in  this 
material  and  considering  all  the  possible  microstructures  present  there. 

2.2.1  Main  result 

With  the  same  setting  as  that  described  for  periodic  homogenization,  we 
may  now  briefly  describe  the  main  result  of  stochastic  homogenization.  The 
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solution  u£  to  the  boundary  value  problem 


-div  (.4  (iu,)  Vu‘)  =  /  in  V,  (5) 

ue  =  0  on  dV, 


converges,  when  e  — >  0,  to  the  solution  uk  of  (4)  where  the  homogenized 
matrix  is  now 


[A*]ij  =  E  ^  jT  (e*  +  Vwei(y,  -))T  ^  (?/»  0  ei  ^  • 


The  corrector  problem  now  reads 

-div[A(y,u)(p  +  Vwp(y,uj))]=0  on 


Vwp  is  stationary,  E  (  /  Vwp(y ,  •)  dy  )  =  0. 

iQ 


(6) 


(7) 


A  striking  difference  between  the  stochastic  setting  and  the  periodic  setting 
can  be  observed  comparing  (2)  and  (7).  In  the  periodic  setting,  the  corrector 
problem  is  posed  on  a  bounded  domain,  namely  the  periodic  cell  Q.  In  sharp 
contrast,  the  corrector  problem  (7)  of  the  random  setting  is  posed  on  the 
whole  space  and  cannot  be  reduced  to  a  problem  posed  on  a  bounded 


domain.  The  reason  is,  condition  E  yj  Vwp(y,  •)  dy J  =  0  in  (7)  is  a  global 

condition.  It  indeed  equivalently  reads,  because  of  the  ergodic  theorem, 

lim  /  Vwp(y,  •)  dy  =  0  for  any  sequence  of  balls  Br  of  radii  It. 

R — >+00  \Br\  Jb 


The  fact  that  the  random  corrector  problem  is  posed  on  the  entire  space 
has  far  reaching  consequences  for  numerical  practice.  Actually,  this  is  proba¬ 
bly  the  main  source  of  all  the  practical  difficulties  of  stochastic  homogeniza¬ 
tion. 


2.2.2  The  direct  numerical  approach 

Practical  approximations  of  the  homogenized  problem  in  random  homoge¬ 
nization  are  not  easily  obtained,  owing  to  the  fact  that  the  corrector  prob¬ 
lem  (7)  is  set  on  the  entire  space.  In  practice,  truncations  have  to  be  con¬ 
sidered,  and  the  actual  homogenized  coefficients  are  only  obtained  in  the 
asymptotic  regime. 
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Although  A*  itself  is  a  deterministic  object,  its  practical  approximation  A*N 
is  random.  It  is  only  in  the  limit  of  infinitely  large  domains  Qn  that  the 
deterministic  value  is  attained  (the  convergence  linpv^oo  A*N(u)  =  A*  has 
been  shown  in  [16,  Theorem  1]). 

At  fixed  A,  the  approximate  homogenized  matrix  A*N  is  random:  a  set  of 
M  independent  realizations  of  the  random  coefficient  A  are  therefore  consid¬ 
ered.  The  corresponding  truncated  problems  (9)  are  solved,  and  an  empirical 
mean  of  the  truncated  coefficients  (8)  is  inferred.  This  empirical  mean  only 
agrees  with  the  theoretical  mean  value  of  the  truncated  coefficient  within  a 
margin  of  error  which  is  given  by  the  central  limit  theorem  (in  terms  of  M). 
For  a  sufficiently  large  truncation  size  N,  this  truncated  value  is  admittedly 
the  exact  value  of  the  coefficient.  The  overall  computation  described  above 
is  thus  very  expensive,  because  each  realization  requires  a  new  solution  to 
the  problem  (9)  of  presumably  large  a  size  since  N  is  taken  large. 


3  Bounds  for  homogenization 

[Work  expanded  in  [1].] 

Given  the  above  computational  workload,  practitioners,  especially  scien¬ 
tists  from  the  applied  communities  (computational  mechanics,  ...),  some¬ 
times  choose  to  avoid  computing  actual  homogenized  equations  and  concen¬ 
trate  on  bounds  on  the  homogenized  matrices  A*.  In  [1],  we  have  carefully 
studied  this  approach,  which  has  some  (strong,  as  will  be  seen  below)  limi¬ 
tations. 
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We  consider  here  the  specific  case  of  composite  materials  consisting  of 
only  two  phases.  We  denote  by  A  and  B  the  associated  matrix  coefficients, 
modelling  the  properties  of  the  phases.  We  also  fix  the  average  volume  frac¬ 
tion  9  of  the  phase  A.  For  simplicity,  we  assume  here  that  9  is  uniform  over 
the  whole  material.  The  problem  is  to  find  all  possible  homogenized  materi¬ 
als,  that  is,  mathematically,  matrices  A* ,  that  can  be  attained  homogenizing 
such  phases  A  and  B  with  the  volume  fraction  9. 

In  this  specific  case,  some  bounds  on  the  homogenized  coefficients  may 
be  established.  Here,  we  present  one  example  of  such  bounds  (actually  the 
most  famous  one).  The  case  we  consider  is  a  scalar  equation  of  the  type  (1) 
with  a  matrix  coefficient  Ae(x)  that  needs  not  be  periodic,  nor  stationary 
ergodic,  and  that  reads 

Ae(x)  =  xe{x)A  +  (1  -  xe(x))B 

where  x£{x)  is  the  caracteristic  function  of  phase  A.  Obtaining  estimates 
on  A *  without  being  in  position  to  explicitly  compute  A *  at  a  reasonable 
computational  price  is  the  whole  interest  of  the  approach  by  “bounds” . 


3.1  The  Hashin-Shtrikman  bounds 

Based  on  the  density  of  the  matrices  obtained  by  periodic  homogenization 
in  the  set  of  matrices  obtained  by  arbitrary  homogenization,  it  is  possible 
to  derive  the  following  Hashin-Shtrikman  bounds  on  A*.  In  the  sequel,  we 
assume  B  >  A. 

Under  the  above  assumptions,  any  homogenized  matrix  A*  satisfies  the 
upper  bound 


A*p  •  p  <  Bp  •  p  +  9  min  [2 p  •  r\  +  (B  —  A)  1r/  ■  r/  +  (1  —  6)h(r])]  (10) 


for  any  p  6  Rd,  where  h(rj)  is  defined  by 


h(r ])  =  min 


1 77  •  fcp 


kezd,k^o  Bk  ■  k 

Similarly,  any  homogenized  matrix  A*  satisfies  the  lower  bound 
A*p  •  p  >  Ap  •  p  +  (1  —  9)  max  [2 p  •  r/  —  (B  —  A)~lrq  •  77  —  9g{rj)\  , 


(11) 
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where  g{rf)  is  defined  by 


g(v) 


\r]  ■  k\2 
max  — — - 

kezd,k^o  Ak  ■  k 


Furthermore,  the  upper  bound  can  always  be  attained:  for  any  p  G  M.d,  there 
exists  a  function  ;y,  Zd-periodic  and  that  generally  depends  on  p,  such  that 
for  the  matrix  A*  obtained  by  periodic  homogenization  of 

A(-J  =  X(-)A+(l-x(-))B, 

£  £  £ 

the  inequality  (10)  becomes  an  equality  (see  e.g.  [30]).  Likewise,  the  lower 
bound  (11)  can  always  be  attained.  We  have  summarized  in  [1,  Section  2.3.2] 
a  proof  of  the  Hashin  Shtrikman  bounds. 


Remark  1  Besides  the  Hashin- Shtrikman  bounds,  many  other  estimates  have 
been  proposed,  such  as  the  dilute  approximation,  the  self-consistent  method  [31] 
and  the  Mori  Tanaka  methods  [28] .  They  are  all  based  on  the  fact  that  the 
problem  of  a  single  inclusion  in  an  infinite  material  (Eshelby  problem)  is  ana¬ 
lytically  solvable  [23].  Similarly  to  the  Hashin- Shtrikman  bounds,  the  spatial 
distribution  of  the  phases  is  not  taken  into  account  in  these  other  bounds. 
The  accuracy  of  these  estimates  and  bounds  strongly  depends  on  the  contrast 
between  A  and  B  and  the  volume  fraction  6  as  shown  on  Figure  1  below. 


3.2  Numerical  illustration 


We  consider  a  two- phase  composite  with  A  and  B.  We  denote  by  a  the 
scalar  conductivity  of  A  (respectively  b  the  conductivity  of  B)  with  a  <  b. 
We  denote  by  d  the  dimension,  and  by  9  the  volume  fraction  of  A. 

We  consider  the  case  of  the  random  checkerboard,  for  which  the  exact 
homogenized  matrix  is  known:  A*  =  a*  Id  =  \/ab  Id.  In  this  simple  case, 
the  different  bounds  and  estimates  have  an  analytical  form:  the  homogenized 
coefficient  a*  is  bounded  from  below  by  the  harmonic  mean  (often  called  the 
Reuss  bound)  and  from  above  by  the  arithmetic  mean  (often  called  the  Voigt 
bound) : 


1 

0/a+(l-0)/b 


<  a*  <  9a  +  (1  -  9)b. 
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These  bounds  are  also  called  Wiener  Bounds  or  Paul  bounds.  In  this  case 
the  Hashin-Shtrikman  bounds  detailed  above  read  (see  e.g.  [25,  page  188]): 


/  d(l-0)(b-a)\ 

v  da  +  0(b-a)  ) 


<  a*  <b 


dd{b  —  a) 
db  +  (1  —  6)  (a 


and  the  Self-Consistent  model  leads  to  an  estimate  Aeg  of  the  effective  con¬ 
ductivity  a*  defined  implicitly  (see  [26])  by 


q  a  ^eff  +  ^  ^eff 

a  +  2Aeff  b  +  2Aeff 


0. 


On  Figure  1  we  plot  these  bounds  and  estimates  for  different  values  of  the 
contrast,  defined  by  b/a ,  for  a  =  1.  Note  that  in  this  case,  by  construction, 
the  volume  fraction  for  any  a  and  b  is  0  =  1/2.  In  Tab.  1,  we  collect  the 
values  of  all  these  bounds  and  estimates,  for  the  particular  case  a  =  1  and 
b=  10. 


Figure  1:  Different  bounds  for  the  checkerboard  test  case. 

A*  Harmonic  HS-  SC  Model  HS+  Arithmetic 
3.16  L81  T38  TOO  4T9  5A0 

Table  1:  Values  of  bounds  and  estimate  for  a  contrast  of  b/a  =  10. 

We  verify  that,  for  the  critical  volume  fraction  0  =  0.5,  even  for  a  contrast 
which  is  not  extremely  large  (a  =  1  and  b  =  10),  the  range  of  homogenized 
matrix  atteignable,  given  by  the  Hashin-Shtrikman  bounds,  is  wide.  In  such 
a  case,  the  spatial  distribution  of  phases,  which  is  not  taken  into  account  on 
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the  bounds,  is  certainly  important.  Note  also  that  a  typical  case  for  real- 
world  composites  is  more  challenging  than  the  case  above,  since  the  contrast 
is  usually  larger  (of  the  order  100  rather  than  10)  and  the  volume  fraction  is 
similar. 

Our  numerical  example  therefore  shows  that,  in  many  cases,  the  Hashin 
Shtrikman  bounds  cannot  provide  accurate  estimates  of  the  homogenized  ma¬ 
trix.  For  a  contrast  of  10,  the  error  between  the  bounds  and  A*  is  larger  than 
25  %.  For  a  contrast  of  100,  the  upper  bound  is  three  times  as  large  as  the 
actual  homogenized  value,  which  is  itself  three  times  as  large  as  the  lower 
bound.  There  is  therefore  a  need  for  developing  efficient  numerical  methods 
that  provide  more  accurate  results. 

4  A  weakly- stochastic  MsFEM  approach 

[Work  expanded  in  [3].] 

In  this  section,  we  show  how  the  Multiscale  Finite  Element  Method  (Ms¬ 
FEM)  can  be  adapted  to  the  stochastic  setting.  We  refer  to  [20]  for  a  review 
on  the  MsFEM  approach.  Let  us  recall  here  that  this  method  is  designed 
to  directly  address  the  original  problem  (namely  (5)  in  the  case  of  interest 
here),  keeping  z  at  its  fixed  value,  rather  than  studying  the  limit  problem 
when  z  — >  0  (as  we  do  in  Section  2,  going  from  (5)  to  (4)).  Another  interest 
of  this  method  is  that  it  does  not  require  any  explicit  formula  for  the  homog¬ 
enized  tensor  (such  as  (2)-(3),  or  (6)-(7)),  which  are  not  always  available. 
More  details  and  comprehensive  numerical  tests  are  published  in  [3].  See 
also  [4], 

4.1  MsFEM  approaches 

For  consistency  and  to  set  our  notation,  we  briefly  review  the  classical,  de¬ 
terministic  setting  for  MsFEM  approaches.  We  next  turn  to  the  stochastic 
setting.  We  consider  problem  (1),  which  we  reproduce  here  for  convenience, 

f  — di v(A£(x)Vu£(x))  =  f(x)  in  V,  (  . 

y  u£  =  0  on  dV, 

where  A£  is  a  symmetric  matrix  satisfying  the  standard  coercivity  and  bound¬ 
edness  conditions.  Note  that  the  approach  is  not  restricted  to  the  periodic 
setting,  so  we  do  not  assume  that  A£(x )  =  A(x/z)  for  a  periodic  matrix  A. 
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As  recalled  above,  we  wish  here  to  keep  e  fixed  at  a  (small)  given  value. 
The  MsFEM  approach  consists  in  performing  a  variational  approximation  of 
(12)  where  the  basis  functions  are  defined  numerically  and  encode  the  fast 
oscillations  present  in  (12).  In  the  sequel  we  will  argue  on  the  variational 
formulation  of  (12): 

Find  u£  E  Hq{T>)  such  that,  Vu  E  Hq(T>),  A£(u£,v )  =  b(v),  (13) 


where 


A£(u,v)  =  y~]  f  A£]{x)<-^-^—dx  and  b(v) 
.  .  Jv  UX{  OXj 


We  introduce  a  classical  PI  discretization  of  the  domain  V,  with  L  nodes, 
and  denote  0°,  i  =  1 , ,L,  the  basis  functions. 


Definition  of  the  MsFEM  basis  functions  Several  definitions  of  the 
basis  functions  have  been  proposed  in  the  literature  (see  e.g.  [24,  21,  22,  14]). 
They  give  rise  to  different  methods.  In  the  following,  we  present  one  of  these 
methods.  We  consider  the  problem 

f  — div(A£(x)V^’K)  =0  in  K 

l  #K  =0°  Ik  ondK. 

Note  the  similarity  between  (14)  and  the  corrector  problem  (2).  Note  also 
that  the  problems  (14),  indexed  by  K,  are  all  independent  from  one  another. 
They  can  hence  be  solved  in  parallel,  using  a  discretization  adapted  to  the 
small  scale  e. 


Macro  scale  problem  We  now  introduce  the  finite  dimensional  space 

Wh  ■■=  span  {4>£,  i  =  1, . . . ,  L}  , 

where  f>\  is  such  that  f f|K  =  (ff  'K  for  all  K,  and  proceed  with  a  standard 
Galerkin  approximation  of  (13)  using  W/, : 

Find  u£h  E  Wh  such  that,  Vw  E  Wh,  A£(u£h, v)  =  b(v).  (15) 

The  dimension  of  VV/t  is  equal  to  L:  the  formulation  (15)  hence  requires 
solving  a  linear  system  with  only  a  limited  number  of  degrees  of  freedom. 
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Numerical  illustration  In  order  to  illustrate  the  MsFEM  approach,  we 
solve  (12)  in  a  one  dimensional  setting  with 

A£(x )  =  5  +  50  sin2  ^ ^  , 

on  the  domain  V  =  (0, 1),  with  s  =  0.025  and  /  =  1000.  We  subdivide  the 
interval  (0, 1)  in  L  =  10  elements.  On  Figure  2,  we  plot  the  MsFEM  basis 
functions  in  a  reference  element  and  the  MsFEM  solution  ueh. 


Figure  2:  Left:  Multiscale  basis  functions  </>e,K  in  the  reference  element. 
Right:  MsFEM  solution  u£h  in  the  domain  (0, 1). 


Natural  adaptation  to  the  stochastic  setting  When  applied  to  the 
stochastic  problem 

f  —div(A£(x,cj)Vu£(x,u;))  =  f(x)  mV,  ,  . 

[  u£  =  0  on  dV, 

where  the  practical  issue  is  to  build  an  estimate  of  the  mean  E {ue(x,  •))  using 
a  Monte-Carlo  simulation  method,  the  natural  adaptation  of  the  MsFEM 
method  is  the  following:  for  each  random  realization  m,  first  construct  a 
MsFEM  basis  and  next  solve  the  macroscale  problem.  This  approach  requires 
a  significantly  large  number  of  computations,  since,  for  each  realization,  a 
new  basis  of  oscillating  functions  is  built,  and  a  problem  at  the  macroscale 
is  solved.  Such  an  approach  has  been  described  and  analyzed  in  e.g.  [18]. 

4.2  A  weakly  stochastic  setting 

As  seen  above,  considering  general  random  materials  lead  to  extremely  ex¬ 
pensive  computations.  The  question  arises  to  know  whether  this  general  ran¬ 
dom  context  is  really  relevant,  and  whether  adequate  modifications  can  lead 
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to  substantial  improvements.  Our  line  of  thought  here  is  based  on  the  fol¬ 
lowing  two-fold  observation:  classical  random  homogenization  is  costly  but 
perhaps,  in  a  number  of  situations,  not  necessary.  Many  materials,  albeit 
not  deterministic,  are  not  totally  random.  Some  of  them  can  be  considered 
as  a  perturbation  of  a  deterministic  material.  The  homogenized  behaviour 
should  expectedly  be  close  to  that  of  the  underlying  deterministic  material 
(and  thus  tractable  from  the  practical  viewpoint),  up  to  an  error  depending 
on  the  amount  of  randomness  present. 

Model  We  introduce  and  study  here  a  specific  model  for  a  randomly  per¬ 
turbed  deterministic  material  (we  refer  to  [12]  for  a  quick  overview  of  this 
setting,  and  some  of  the  associated  numerical  techniques,  and  to  [1]  for  a 
more  comprehensive  review  of  our  contributions  along  these  ideas).  We  are 
interested  in  the  following  elliptic  problem 

f  —  div  (A*  (a;,  u>)  Vu®)  =  f(x)  in  V  C  M.d, 

[  =  0  on  dV, 

that  is  (16)  with 

A£(x,lj)  =  A£v(x,uj)  =  A£0{x )  +  r/AKx,!^), 

where  i]  G  M  is  a  small  parameter,  Aq  is  a  deterministic  matrix  uniformly 
elliptic  and  bounded,  and  Af  (x,  u>)  is  a  bounded  matrix.  The  matrix  A®  is 
hence  a  perturbation  of  the  deterministic  matrix  Af,. 

Remark  2  The  above  setting  is  of  course  one  possible  setting  where  the  the¬ 
ory  may  be  developed.  Other  forms  of  random  perturbations  of  deterministic 
( possibly  periodic)  problems  could  also  be  addressed.  See  e.g.  [5,  6,  1,  8]  and 
the  review  article  [1]. 

In  the  case  (18),  a  MsFEM  method  alternative  to  the  one  presented  in 
Section  4.1  can  be  proposed.  The  idea  is  to  compute  the  MsFEM  basis 
functions  only  once ,  with  the  deterministic  part  of  the  matrix  A);  and  next 
to  perform  Monte-Carlo  realizations  only  for  the  macro  scale  problems.  We 
refer  to  [3]  for  all  the  details. 

As  above,  we  hence  first  solve  (14),  with  As  =  Af ,  and  build  the  finite 
dimensional  space 


(17) 

(18) 


Wh  :=  span {$,  i=l,...,L}. 
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We  next  proceed  with  a  standard  Galerkin  approximation  of  (17)  using  W/p 
for  each  m  E  {1, . . . ,  M},  we  consider  a  realization  and  compute 

um(-,Lj)  G  Wh,  such  that 

_  r  f)v  r 

Vv  G  Wh,  J  (Avm)i3(xiu)-^{x,Lj)  —  (x)dx  =  J  f(x)v(x)dx.  (19) 

Since  the  MsFEM  basis  functions  are  only  computed  once,  a  large  computa¬ 
tional  gain  is  expected,  and  this  is  indeed  the  case. 

Numerical  studies  We  now  estimate  the  performance  of  the  approach.  To 
this  aim,  we  compare  the  solution  of  the  above  approach  with  the  solution  of 
the  direct  approach  (of  Section  4.1)  and,  for  reference,  the  solution  to  (17) 
obtained  using  a  finite  element  method  with  a  mesh  size  adapted  to  the  small 
scale  e.  Our  estimators  are  built  as  follows: 

e(u!,  u2)  =  E  n^iJ1  )  ,  (20) 

where  N  is  the  norm  used,  u  i  and  u2  are  the  solutions  obtained  with  any  two 
different  methods.  The  expectation  is  in  turn  computed  using  a  Monte-Carlo 
method.  Considering  M  realizations  {Xm{u)}l<m<M  of  a  random  variable 

(here  X(u)  =  -  —  ’  \ - ttp — ^— ),  we  compute  the  empirical  mean  uM 

IMuWlIv 

and  the  empirical  standard  deviation  gm  as 
Vm(X) 

As  a  classical  consequence  of  the  Central  Limit  Theorem,  it  is  commonly 
admitted  that  E(X)  satisfies 

E(X)-iim(X)\  <  1.96 

We  consider  the  following  numerical  example.  Let  (X^)^  ^eZ2  denote  a 
sequence  of  independent,  identically  distributed  scalar  random  variables  uni¬ 
formly  distributed  over  the  interval  [0, 1].  We  define  the  random  conductivity 


m= 1 

i  m 

m= 1 
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matrix  as 


A£v(x,y,u) 


2  +  Psin(27ra;/£) 
2  +  P  sm(2ny  /  e) 


2 -f  sm(2ny / s)  \ 
2  ±  P  sm(27vx  /  e) ) 


(1  Pr]Xkti(uj))ld2, 


with  the  parameters  P  =  1.8  and  e  =  0.025.  Then  we  compute  uref  solu¬ 
tion  to  (17)  on  the  domain  V  =  (0,  l)2,  with  /  =  1.  Let  um  and  us  be 
its  approximation  by  the  general  MsFEM  approach  (of  Section  4.1)  and  the 
weakly-stochastic  MsFEM  approach  described  above,  respectively.  The  nu¬ 
merical  parameters  for  the  computation  are  determined  using  an  empirical 
study  of  convergence.  We  used  for  the  reference  solution  a  fine  mesh  of  size 
hf  =  e/40.  The  MsFEM  basis  functions  are  computed  in  each  element  K 
using  a  mesh  of  size  Km  =  e/80.  The  coarse  mesh  size  is  H  =  1/30.  We 
consider  M  —  30  realizations. 

We  report  in  Tables  2  and  3  the  estimator  (20),  along  with  its  confidence 
interval,  for  the  norms  Hl{V)  and  L2(V),  respectively. 


V 

^  'LL  ref) 

ciflls-)  U ref  ) 

e(us,  um) 

1 

0.1 

0.01 

8.12  ±0.19 
7.17  ±0.02 
7.15  ±0.002 

17.37  ±0.78 
7.62  ±0.07 
7.28  ±0.007 

15.51  ±0.87 
2.56  ±0.10 
1.39  ±0.002 

Table  2:  Hl(T>)  error  (in 

%). 

V 

^  'U'ref) 

^ ref  ) 

e(us,uM) 

1 

0.1 

0.01 

0.56  ±0.08 

0.54  ±0.01 
0.53  ±0.001 

1.69  ±0.49 

0.57  ±0.06 
0.62  ±0.005 

1.47  ±0.50 
0.20  ±0.07 
0.11  ±0.003 

Table  3:  L2(V)  error  (in  %). 


We  observe  that,  when  rj  is  small  (here,  r/  <  0.1),  the  alternative  approach 
provides  a  function  us  that  is  an  approximation  of  uref  as  accurate  as  %. 
Recall  that,  since  the  MsFEM  basis  has  only  been  computed  once,  the  cost 
for  obtaining  an  empirical  approximation  of  is  much  smaller  than  that 
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for  obtaining  the  corresponding  empirical  estimator  of  E(«j^).  This  demon¬ 
strates  the  efficiency  of  the  approach.  As  expected,  when  rj  is  not  small 
(say  rj  ~  1),  the  accuracy  of  the  solution  ug  computed  with  the  alternative 
approach  proposed  here  decreases. 

Elements  of  proof  In  [3],  we  have  analyzed  the  method  introduced  here 
in  the  one-dimensional  setting  (see  also  [4]).  For  the  sake  of  analysis,  we 

assume  that  the  highly  oscillating  coefficient  reads  a£v(x,uj )  =  av  .uj, 

where  av  satisfies  the  standard  assumption  of  stochastic  homogenization  (see 
Section  2.2).  The  problem  (17)  now  reads 


We  assume  that  the  randomness  is  small,  in  the  sense  (see  (18))  that 

av(x,  oS)  =  aper(x)  +  r/  ai(x,  u),  (22) 

where  aper  is  a  deterministic,  periodic  function  and  r/  is  a  small  deterministic 
parameter. 

In  [3],  we  have  bounded  from  above  the  difference  between  u*,  the  so¬ 
lution  to  the  homogenized  equation  (4),  and  the  weakly-stochastic  MsFEM 
solution,  in  the  following  sense.  For  a  given  realization  of  the  random  co¬ 
efficient  ari(x,u),  let  u(-,u>)  be  the  weakly-stochastic  MsFEM  solution,  that 
solves  (19).  By  construction,  this  solution  is  a  linear  combination  of  the 
highly  oscillating  basis  functions: 

L 

u(x,uj)  =  ^2Ui(uj)<l>i(x). 

2=1 

Let  uw_msfem(f  u;)  be  the  associated  representation  in  terms  of  standard  PI 
basis  functions: 

L 

Av— MsFEM  OF  u)  =  Uj(uj)(f)^(x). 

2=1 

We  have  the  following  estimate: 

K,  -  'V  MsFF.M ( '•  -x')  | // 1  (0. i )  <  C  (h  +  |  +r]a£h(uj)  +  r]2C(r])^  (23) 
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where  C  is  a  deterministic  constant,  independent  of  h,  e  and  //,  and  C{rj)  is  a 
deterministic  function,  bounded  when  77  — ^  0.  In  the  above  bound,  07; (a;)  is 
a  random  number,  independent  of  77,  that  depends  on  s,  h  and  the  random 
term  ai(x,u>)  in  (22). 

Let  us  comment  on  (23).  Assume  that  77  =  0,  i.e.  the  problem  (21)  is  a 
periodic  problem.  Then  our  method  is  identical  to  the  standard  deterministic 
MsFEM  method,  and  we  recover  from  (23)  the  classical  bound  known  in  that 
case,  namely 

\\K)  “  ^,MsFEM||iJ1(0,l)  <  C  (h  +  —  j  . 

Assume  now  that  a\  is  deterministic.  Then  our  method  is  not  exactly  the 
MsFEM  method,  since  we  do  not  take  into  account  a  \  to  build  the  highly 
oscillating  basis  functions.  In  that  case,  cr^(u;)  turns  out  to  vanish,  and  we 
infer  from  (23)  that 

\\u*  —  Av— MsFEM 1 1 h1  (0,1)  <  C  (h  +  —  +  rj2C{r])^  . 

We  hence  observe  that,  provided  the  term  neglected  to  build  the  basis  func¬ 
tions  is  small  (namely  77  <C  1),  we  obtain  a  similar  accuracy  as  with  the 
standard  MsFEM  method. 

A  similar  conclusion  holds  in  the  general  case  (22).  Note  also  that  the 
bound  (23)  is  valid  for  any  realization  u  of  the  randomness.  It  is  therefore 
a  more  precise  result  than  a  bound  on  the  expectation  of  the  error,  where 
all  random  realizations  are  averaged.  For  instance,  the  bound  (23)  allows  to 
understand  what  is  the  probability  distribution  of  the  error. 


5  Conclusions  and  plan  for  the  following  years 

In  this  report,  we  have  first  reviewed  an  approach  to  obtain  bounds  (here,  the 
Hashin-Shtrikman  bounds)  on  the  homogenized  matrix.  This  approach  only 
involves  computations  of  moderate  difficulty.  However,  we  have  outlined  the 
strong  limitations  of  such  a  strategy.  In  some  cases,  the  difference  between 
the  lower  and  upper  bounds  is  indeed  very  large.  The  obtained  estimations 
are  then  inaccurate.  This  motivates  the  development  of  efficient  numerical 
methods  that  provide  more  accurate  results. 

To  this  aim,  we  have  focused  on  weakly  stochastic  materials,  for  which  we 
successfully  adapted  the  well-known  Multiscale  Finite  Element  Method  (Ms¬ 
FEM).  We  have  proposed  a  method  with  a  much  smaller  computational  cost 
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than  the  original  MsFEM  in  the  stochastic  setting.  Provided  the  stochastic 
perturbation  is  indeed  small,  the  method  we  propose  is  as  accurate  as  the 
original  one. 

We  summarize  now  the  directions  of  research  we  wish  to  pursue  during 
the  following  years. 

A  variant  of  classical  random  homogenization  In  the  short  term,  our 
aim  is  to  study  a  particular  setting  for  stochastic  homogenization,  which  is 
not  the  classical  setting  described  in  Section  2.2  (where  the  random  coefficient 
A  in  (5)  is  stationary).  The  setting  we  wish  to  study  is  the  case  when  the 
random  coefficient  is  the  composition  of  a  standard  deterministic  and  periodic 
function  Aper  with  a  stochastic  diffeomorphism: 

A£(x,oj)  =  Aper  <f>-1  (24) 

where,  for  any  random  realization  u>,  the  application  x  i— >  u>)  is  a  diffeo¬ 

morphism.  Formally,  such  a  setting  models  a  microstructure  that  is  periodic, 
in  a  given  reference  configuration.  The  latter  is  only  known  up  to  a  certain 
randomness.  Materials  we  have  in  mind  are  ideally  periodic  materials,  where 
some  random  deformation  has  been  introduced,  for  instance  during  the  manu¬ 
facturing  process.  Othewise  stated,  these  are  periodic  materials  seen  through 
random  glasses!  This  setting  has  been  initially  introduced  in  [9],  where  the 
homogenized  problem  is  identified. 

An  interesting  question  in  that  context  is  that  of  numerical  discretization. 
In  the  classical  context,  a  standard  procedure  is  to  solve  the  corrector  problem 
on  a  truncated  domain  (see  (8)  and  (9)).  The  convergence  of  the  procedure 
is  given  by  [16,  Theorem  1].  We  currently  work  on  a  similar  analysis  in  the 
context  of  (24)  (see  [2]). 

A  more  theoretical  question  is  to  precisely  understand  the  behaviour  of 
u£(x,u)  —  u*(x)  when  £  goes  to  0,  where  u£  is  the  solution  of  the  highly 
oscillating  problem  and  u*  the  solution  of  the  homogenized  problem.  In 
the  classical  setting,  and  in  the  one-dimensional  case,  the  convergence  of 
£~1^2(ue(x,u)  —  u*(x ))  to  a  Gaussian  process  has  been  shown  in  [17].  This 
question,  in  the  context  of  (24),  is  addressed  in  [2], 

The  setting  (24)  is  in  general  not  a  weakly  stochastic  setting,  as  the 
amount  of  randomness  present  in  $  may  be  large.  Yet,  in  the  case  when  the 
diffeomorphism  $  is  close  to  the  identity,  namely 

$(a:,  to)  =  x  +  7]^(x,  uj)  +  0(rj2)  (25) 
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for  a  small  deterministic  parameter  77,  the  amount  of  randomness  turns  out 
to  be  small.  This  case  is  thus  another  instance  of  randomly  perturbed  deter¬ 
ministic  materials  (recall  Section  4.2,  where  we  introduced  another  weakly 
stochastic  setting,  and  Remark  2,  where  we  pointed  out  other  weakly  stochas¬ 
tic  settings).  The  case  (24)-(25)  has  been  studied  in  [10,  11]. 

A  Fast  Fourier  Transform  approach  In  the  course  of  our  investigations, 
we  have  identified  the  following  tracks  of  research,  which  are  closely  related 
to  the  research  directions  of  the  contract. 

First,  in  the  periodic  homogenization  setting  recalled  in  Section  2.1,  a 
method  based  on  Fast  Fourier  Transform  has  been  proposed  in  [29,  27].  The 
idea  is  as  follows.  Let  A0  be  a  constant  symmetric  positive  matrix.  The 
corrector  problem  (2)  is  equivalent  to 

(  -div  (A0  (p  +  Vwp))  =  div((Aper(y)  -  A0)  (p  +  \7wp))  in  Rd, 

1  wp  is  Zd-periodic. 

The  idea  of  [29,  27]  is  to  solve  this  problem  iteratively.  Knowing  the  iterate 

at  iteration  k,  the  next  iterate  wp+1  is  defined  as  the  unique  solution  to 

f  -div  (A0  (p  +  Vu)p+1))  =  div  (( Aper(y )  -  A0)  (p  +  Vw£))  in  Rd, 

1  Wp+1  is  Zd-periodic. 

As  A0  is  a  tensor  independent  of  y  (in  contrast  to  Aper(y)),  the  above  problem 
can  be  solved  very  efficiently  using  a  Fourier  transform.  Hence,  rather  than 
solving  (2),  we  are  left  with  solving  many  times  a  simpler  problem. 

Our  aim  is  to  compare  this  iterative  method  with  the  standard  method, 
in  term  of  efficiency.  The  choice  of  A$  is  most  probably  of  paramount  impor¬ 
tance,  since  the  convergence  rate  (and  also  the  fact  that  the  iterations  in  k 
converge  or  not)  depends  on  it.  We  have  already  run  some  preliminary  tests 
with  this  method,  but  definite  conclusions  are  yet  to  be  obtained. 

Second,  in  the  context  of  stochastic  homogenization,  approaches  using 
some  decomposition  of  the  random  matrix  A(x,  u>)  in  (5)  would  be  worthwhile 
to  investigate. 
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